Analytical model for flux saturation in sediment transport 
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The transport of sediment by a fluid along the surface is responsible for dune formation, dust en- 
trainment and for a rich diversity of patterns on the bottom of oceans, rivers, and planetary surfaces. 
Most previous models of sediment transport have focused on the equilibrium (or saturated) particle 
flux. However, the morphodynamics of sediment landscapes emerging due to surface transport of 
sediment is controlled by situations out-of-equilibrium. In particular, it is controlled by the satura- 
tion length characterizing the distance it takes for the particle flux to reach a new equilibrium after 
a change in flow conditions. The saturation of mass density of particles entrained into transport and 
the relaxation of particle and fluid velocities constitute the main relevant relaxation mechanisms 
leading to saturation of the sediment flux. Here we present a theoretical model for sediment trans- 
port which, for the first time, accounts for both these relaxation mechanisms and for the different 
types of sediment entrainment prevailing under different environmental conditions. Our analytical 
treatment allows us to derive a closed expression for the saturation length of sediment flux, which 
is general and can thus be applied under different physical conditions. 

PACS numbers: 45.70.-n, 47.55.Kf, 92.40.Gc 



I. INTRODUCTION 

When a sediment bed is exposed to a fluid flow, parti- 
cles can be entrained and transported by different mech- 
anisms. The transport regime depends primarily on 
the inertial characteristics of the particles and the fluid. 
Sufficiently light particles are transported as suspended 
load in which their weight is supported by the turbu- 
lence of the fluid. In contrast, particles which are suf- 
ficiently heavy are transported along the surface [1, 2]. 
This type of transport incorporates two main transport 
modes, namely: saltation, which consists of sediment 
grains jumping downstream close to the ground at nearly 
ballistic trajectories, and creep, which consists of parti- 
cles rolling and sliding along the sediment bed. Sedi- 
ment transport along the surface is responsible for a wide 
range of geophysical phenomena, including surface ero- 
sion, dust aerosol emission, and the formation and mi- 
gration of dunes [1-6]. Therefore, the quantitative un- 
derstanding of sediment transport may improve our un- 
derstanding of river beds evolution [2], the emission of 
atmospheric dust [4, 6] and the dynamics of planetary 
sand landscapes [3, 6, 7]. 

Once sediment transport begins, the fluid loses mo- 
mentum to accelerate the particles as a consequence of 
Newton's second law (the transport-flow feedback, e.g. 
[8, 9]). Therefore, the sediment flux, Q, which is the av- 
erage momentum of grains transported per unit soil area, 
is limited by an equilibrium value, the saturated flux, Q s . 
Although previous studies focused on this equilibrium 
flux (e.g. [10-15]), the dynamics of sediment landscapes 
is controlled by situations out-of-equilibrium. In particu- 
lar, the sediment flux needs a spatial lag — the so-called 



saturation length, L s — to adapt to a change in flow 
conditions [16-19]. This saturation length introduces the 
main relevant length-scale in the dynamics of sediment 
landscapes under water and on the surface of planetary 
bodies. For instance, the saturation length controls the 
minimal size of crescent-shaped (barchan) dunes moving 
on top of bedrock, as well as the wavelength of the small- 
est dunes (the so-called "elementary dunes") emerging 
on top of a sediment bed [17, 19]. Although important 
insights were gained recently from experimental studies 
[17, 20, 21], the physics behind the saturation length, 
and thus the dependence of L s on flow and sediment at- 
tributes, is still insufficiently understood. 

One of the most important deficiencies in our under- 
standing of the dependence of L s on flow and sediment 
conditions is that it remains uncertain which mechanisms 
are most important in determining the saturation of the 
sediment mass flux. On the one hand, it has been sug- 
gested that the acceleration of transported particles due 
to fluid drag is the dominant relaxation mechanism [17- 
20]. This model neglects the entrainment of sediment 
bed particles due to fluid lift, as well as the entrainment 
of sediment bed particles and the deceleration of trans- 
ported particles due to collisions of transported parti- 
cles with the sediment bed (grain-bed collisions) . On the 
other hand, the entrainment of sediment bed particles by 
fluid lift and grain-bed collisions has also been proposed 
to be the dominant relaxation mechanism [16, 22]. How- 
ever, this model neglects momentum changes of trans- 
ported particles, which is exactly the opposite situation 
of the models in Refs. [17-20]. Moreover, to our knowl- 
edge, all previous models neglected a further relaxation 
mechanism of the sediment flux, namely, the relaxation 
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of the fluid speed in the transport layer (U) due to the 
saturation of the transport- flow feedback [23] . 

To address this situation and develop an accurate ex- 
pression for L s that can be used in future studies, this 
paper presents a model for flux saturation in sediment 
transport which, for the first time, accounts for all afore- 
mentioned mechanisms for the saturation of sediment 
flux. In particular, our theoretical model accounts for 
the coupling between the entrainment of sediment bed 
particles due to fluid lift and grain-bed collisions, the ac- 
celeration and deceleration of transported particles due 
to fluid forces and grain-bed collisions, and the saturation 
of U due to the saturation of the transport-flow feedback. 
Our analytical model allows us to derive a closed expres- 
sion for L s which can be applied to different physical 
environments. Our model suggests that grain-bed colli- 
sions, which have been neglected in all previous studies, 
have an important influence on the saturation length, L s . 
Moreover, our model suggests that the relaxation of U 
plays an important role for sediment transport in dilute 
fluids (aeolian transport), whereas it plays a negligible 
role for sediment transport in dense fluids (subaqueous 
transport). 

This paper is organized as follows. Sections II and III 
discuss the analytical treatment of flux saturation. In 
the former Section, we derive the mass and momentum 
conservation equations for the layer of sediments in trans- 
port, as well as the differential equation of the sediment 
flux in terms of the mass density and average velocity 
of the transported particles. These equations allow us 
to obtain a mathematical expression for the saturation 
length of sediment transport, which is presented in Sec- 
tion III. This Section also discusses how to determine 
the quantities appearing in the saturation length equa- 
tion, which encode the attributes of sediment and flow, 
as well as the characteristics of sediment entrainment and 
particle-fluid interactions. In Section IV we use our the- 
oretical expression to perform a study of the saturation 
length as a function of the relevant physical quantities 
controlling saturation of sediment flux. Conclusions are 
presented in Section V. 



inverse first-order Taylor coefficient of T(Q), 



L s = - 



ary 1 



(2) 



In this Section, we derive the equations that describe 
the downstream evolution of the sediment mass flux, Q, 
towards its equilibrium value, Q Sl in sediment transport 
under turbulent boundary layer flow. 

The mass flux Q is defined as Q = MV, where M is 
the average transported mass per unit soil area and V is 
the average particle velocity. Therefore, the saturation of 
Q is dictated by the mechanisms governing the relaxation 
of M and V towards their saturated values, M s and V s , 
respectively. The quantitative description of the satura- 
tion processes of M and V requires incorporation of all 
relevant forces acting on the sediment particles in trans- 
port, namely drag, gravity, buoyancy, collision forces be- 
tween particles in transport ("mid-fluid collisions") and 
friction due to collisions between particles and the bed. 
Indeed, Moraga ct al. [24] found experimentally that lift 
forces due to shear flow acting on a particle surrounded 
by fluid — which have often been assumed to be signif- 
icant during transport (e.g. [2]) — are approximately 
an order of magnitude smaller than the drag force and 
can be, thus, neglected in our calculations. On the other 
hand, the so-called added mass force exerted by acceler- 
ated or decelerated particles to dislodge the fluid as they 
move through it leads to enhanced inertia of the parti- 
cles in transport. This added mass effect plays a relevant 
role for the motion of the particles [25] , and thus we also 
take it into account. Our analytical treatment applies to 
situations where the fluid velocity is not too high such 
that only transport through saltation or creep (the main 
transport modes of particles along the surface [1, 6]) is 
considered. Transport through suspension or in dense 
transport regimes, such as sheet flow [26], are, thus, not 
considered. 

In Section II A we first present the definitions and no- 
tations used in our study. Afterwards in Section II B we 
derive the equations describing the kinematics of trans- 
ported grains, from which we obtain the saturation equa- 
tions, presented in Section II C. 



II. FLUX SATURATION IN SEDIMENT 
TRANSPORT 



Definitions and Notations 



The downstream evolution of the sediment flux, Q, to- 
wards its equilibrium value, Q s , can be described by the 
following equation [18], 



r(Q) 



dQ _ Qs-Q 
dx L, 



(1) 



which is valid in the regime where Q is close to saturation 
(|1 — Q/Q s \ <C 1)- The length-scale L s , the saturation 
length, characterizes the response of the sediment flux 
due to a change in flow conditions around the equilib- 
rium. Since r(Q s ) = 0, L s can be written as the negative 



We use a three-dimensional coordinate system (x, y, z), 
where x denotes the direction of fluid motion, y is the 
lateral direction and z is the vertical direction. The top 
of the sediment bed, which corresponds to the height 
at which the local particle concentration equals approxi- 
mately 50% of the particle concentration deep within the 
bed [27], is located at the vertical position z = h D (x,y). 
Particles at z > h a are transported by fluid forces, 
whereas transport of particles at height z below h D is 
neglected. 

Since the time-scale of the relaxation of the sediment 
flux due to changes in the flow is typically much smaller 
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than the time-scale of the evolution of bedforms (dunes 
and ripples) [19], we can adopt the approximation that 
transport over the sediment landscape is in the steady- 
state, i.e. d/dt = 0, where t denotes time. The positions 
x = (x, y, z) and velocities v = (v x ,v y ,v z ) of individual 
transported grains (z > h ) are described using the phase 
space density, 



dN = P(x,v)d 3 xd 3 v, 



(3) 



which gives the number of grains dN moving in the in- 
finitesimal phase space volume d 3 a;d 3 w. P(x, v) is a func- 
tion which maps the space M 6 : z > h (x,y) to K, and 
which must vanish if z — > oo or |v| — > oo. In particular, 
there must be a finite cut-off height z cu t < oo and a finite 
cut-off velocity u cut < oo above which P vanishes, 



z > z C ut : P = 0, 
|v| > v cut : P = 0. 



(4) 
(5) 



This is because P is a continuum approximation of an 
actually discrete system with a finite number of particles 
moving above a fixed area of the sediment bed, and a set 
with a finite number of elements has a largest element; 
in other words, P as a function of z and v has a compact 
support. Note that the cut-off height is typically much 
smaller than the height of the turbulent boundary layer, 
8b, 



Z C ut < 6 b - 



(6) 



The mass per unit phase space volume is computed as 
mP, where m = p p ^d 3 is the mass of each particle with 
p p the particle density and d the mean grain diameter. 
Consequently, the local mass of transported particles per 
unit volume, p, is given by the equation, 



PdV 



(7) 



In the calculations which follow in the next Sections, 
average values of a quantity / are computed according 
to the following definitions, 



/ Pfd 3 v 



m 
9 



J Pd 3 v 



- / Pfd'v 



(fh = 



J Pfd 3 v 
/ Pd 3 v 



= - I Pfd\ 



</h = 



/ Pfd 3 v 
/ Pd 3 v 



J P fd2 



-7 h a 



M I Pfd2 > 



(12) 



(13) 



= - / Pfd 3 v; (14) 
Pi 



(15) 



where we used Eqs. (7-11). By using Eqs. (12) and (15), 
we obtain the equations for the total mass flux Q and the 
average particle velocity V, which read, respectively, 



Q = J p(v x )dz = M(v x ), 

h 



(16) 
(17) 



while the total transported mass per unit soil area, M, 
is calculated from the expression, 



B. Kinematics of the particles transported by the 

fluid 



M = J pdz. 



(8) 



Analogously, the local mass of transported upward- 
moving (pf) and downward-moving (pi) particles per 
unit volume are given by, 



(9) 
(10) 



p^ = m I Pd v, 



Pi = m / Pd v, 



P = Pt + Pl, (11) 
where Ki = R 3 : v z > and R 3 = R 3 : v z < 0. 



In this Section, the kinematics of the transported 
grains arc described using the notations and definitions 
introduced in the last Section. First, the Boltzmann 
equation for the flux of particles in transport is derived, 
whereupon we derive the mass and momentum conserva- 
tion equations of the particles transported by the fluid. 



1. Boltzmann equation for sediment flux 

The kinematics of the transported particles are de- 
scribed in terms of the Boltzmann equation, which writes, 
in the most general form [28, 29], 
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where V v = (d/dv x ,d/dv y ,d/dv z ). This equation states 
that the phase space density P of transported particles, 
as defined by Eq. (3), changes due to non-conservative 
forces [28] and collisions between particles, which occur 
with a rate given by the function Ki(P). The mean free 
path of a transported particle, i.e. the average down- 
stream distance a particle travels without colliding with 
any transported particle, is typically much larger than 
the grain diameter — as stated previously, dense trans- 
port regimes as sheet flow are not considered in our treat- 
ment of flux saturation. By taking this difference in 
length-scales into account, the particles can be consid- 
ered as points, and it can be shown that [30], 



Ki(P)d s v = 0, 
vK l (P)d 3 v = 0. 



(19) 
(20) 



These equations express that the total number and mo- 
mentum of colliding particles is conserved. It is impor- 
tant to emphasize that the Boltzmann equation is usu- 
ally applied for systems under the influence of conserva- 
tive external forces (e.g. [29]). However, in our system, a 
transported grain with position x and velocity v traveling 
at a height z > h a is subjected to a non-conservative force 
F(x, v), which incorporates gravity, buoyancy, drag and 
added mass, while mid-fluid collisions are incorporated 
in the function Ki(P). Since for saltation and creep, 
the effect of time-fluctuations of the fluid velocity profile 
on the particle motion is small [31], time-fluctuations of 
F(x, v) can be neglected. In this manner, dP/dt can be 
written as, 



dP F 

— = v-VP+--V v P, 

at m 



(21) 



where V = (d/dx,d/dy,d/dz) while we used d/dt = 
0. Hence, the Boltzmann equation (Eq. (18)) for the 
phase space density of the particles in transport above 
h„ writes, 



v • VP + V v 



(22) 



2. Mass and momentum, conservation equations 

In order to derive the mass and momentum conserva- 
tion equations for the flux of particles transported by the 



fluid, we use the following results, 



V v • (PF)d 3 u = / (PF) ■ ndS = 



(23) 



R3 



vV v • (PF)d 3 w = / (PvF) • ndS 



(V v v)-(PF)d 3 w = -i-(F). (24) 
m 



R3 



In order to obtain the equations above, we performed 
partial integration and used the Stokes theorem, and we 
further used Eq. (12) and V v v = I (the identity tensor). 
Furthermore, we also took into account that, since P has 
a compact support (see Eq. (5)), the integrals over the 
boundary <9R 3 of the velocity space vanish. 

Mass conservation equation - - we then integrate 
Eq. (22) over the whole velocity space (J" R3 -d 3 i>) in or- 
der to obtain the mass conservation equation. By using 
Eqs. (12) and (23), and by noting that the integral of 
the right-hand-side of Eq. (22) over the velocity space is 
equal to zero due to Eq. (19), we obtain, 



V-(p(v»=0. 



(25) 



Momentum conservation equation — next, we multiply 
Eq. (22) with v and integrate the resulting equation over 
the whole velocity space (J K3 -vd 3 u) in order to obtain 
the momentum conservation equation. Hereby we use 
Eqs. (12) and (24), and the fact that the integral of the 
right-hand-side of Eq. (22) multiplied with v over the 
velocity space vanishes due to Eq. (20). In this manner, 
we obtain, 



V • (p(vv)) 



(F). 



(26) 



Equations (25) and (26) are the mass and momentum 
conservation equations, respectively, of the flux of parti- 
cles transported by the fluid, i.e. those particles traveling 
at height z > h a . 



C. Differential equations of flux saturation 

The results of the last Section can be now used in or- 
der to derive the saturation equations for the average 
transported mass per unit soil area (M) and the average 
particle velocity (V), used to define the sediment flux, 
Q = MV. To do so, we first integrate Eqs. (25) and (26) 
over the height. This calculation is the subject of Section 
II C 1. Thereafter, in Section II C 2 we rewrite each term 
of the resulting momentum balance equation in terms of 
M and V. We then present the mass and momentum 
balance equations in their final form in Section II C 3. 



5 



1. Height-integrated conservation equations 

Since our description relates to the saturation of the 
mass flux Q due to changes in the downstream direc- 
tion, we can integrate Eq. (25) and the horizontal com- 
ponent of Eq. (26) over the height (J, -dz) and consider 
a laterally invariant sediment bed (d/d y = 0). By using 
Eqs. (15)-(17), and by further taking into account that 
p(z) also has a compact support like P does (cf. Eq. (4)), 
the height-integration of the horizontal components of 
Eqs. (25) and (26) yields, 

d -(MV) = ( P (v z ))(ho), (27) 



dx 



dx 
(c v MV 2 ) 



M 



m 



(F x ) + {p{v x v z ))(h ), (28) 



where, 



Cy 



V 2 



(29) 



is a correlation factor. Equation (27) means that the 
change of the mass flux Q = MV with x equals the 
amount of material which enters or leaves the system 
at the top of the sediment bed {{p(v z ))(h )) due to ero- 
sion and deposition. Further, Eq. (28) means that the 
change of the energy density ^M(v 2 ) = ^c v MV 2 with x 
is governed by the inc rease of particle momentum due to 
fluid acceleration (^(F x )) and the loss of particle mo- 
mentum due to friction at the top of the sediment bed 
{{p(v x v z )){h )). 

The correlation factor — we note that c v is invari- 
ant under the scaling transformation P — > c\P, where 
ci is a constant factor. For transport in the equilib- 
rium (d I 'dx — 0), c v is nearly constant with the fluid 
shear velocity, u*. In fact, experiments and numerical 
simulations on sediment transport in the aeolian regime 
[6, 9, 27, 32, 33] show that both p(z)/M and (v x )(z) 
- and, as a consequence of Eq. (29), c v — are nearly 
independent of it*. On the other hand, numerical simu- 
lations on sediment transport in the subaqueous regime 
[27] show that, in this regime, p(z) and (v x )(z) change 
nearly scale-invariantly with u*, i.e. p(z) — > C2p(z), 
( v x){z) cz{vx){z), where ci and C3 are constant fac- 
tors. Therefore, P also changes nearly scale-invariantly 
with it*, such that c v is nearly constant with also for 
sediment transport in the subaqueous regime. By con- 
sidering that c v is nearly independent of u* , it seems rea- 
sonable that changes of c v with x during the saturation 
process of the sediment flux close to the equilibrium can 
be regarded as negligible compared to the corresponding 
changes of M or V with x. In this manner, we can con- 
sider the value of c v associated with sediment transport 
in the equilibrium, independent of the downstream po- 
sition and of the fluid shear velocity. This leads to the 
following approximation for c„ , 



(30) 



where (v 2 ) s is the equilibrium value of (v x ). The equilib- 
rium value of c v can be determined from experiments as 
we will discuss in Section III C 1 and III D 1. 



2. Momentum balance equation in terms of M and V 

Now we express both terms on the right-hand-side of 
the momentum conservation equation, i.e. Eq. (28), as 
functions of M and V in order to obtain a differential 
equation describing the saturation of M and V. 

The first term on the right-hand-side of Eq. (28) can 
be written as [9, 16], 



M 



(Fx 



3M 
Ac a sd 



C d (V r ) ■ V, 



2 

r ' 



(31) 



where s = p p /pf is the ratio between sediment and fluid 
density; V r is defined as, 



V r = U — V, 



(32) 



which is the difference between the average fluid velocity 
(U = u) and the average horizontal particle velocity (V), 
where u(z) is the fluid velocity profile; Cd is the drag 
coefficient, which is a function of V r , and c Q accounts for 
the added mass force through, 



1 



1 

2s' 



(33) 



v 2 



The added mass force arises when the particle is acceler- 
ated relative to the surrounding fluid, because the fluid 
layer immediately surrounding the particle will also be 
accelerated. As denoted by Eq. (33), this "added mass" 
of the fluid layer amounts to approximately one half the 
weight of the fluid displaced by the particle [2]. While 
the added mass correction is significant for transport in 
a dense medium such as water [25], it is negligibly small 
for sediment transport in the aeolian regime since c a ~ 1 
for large s. Thus, this correction is usually disregarded in 
studies of aeolian sediment transport (e.g. Refs. [9, 16]). 
We note that Eq. (31) is not valid for dense transport 
regimes like sheet flow, in which the drag coefficient dis- 
plays a strong dependence on the concentration profile of 
transported particles [34]. In this manner, Eq. (31) can 
be used in the present study because our analytical treat- 
ment considers the two main modes of transport, namely 
saltation and creep. 

The second term on the right-hand-side of Eq. (28) is 
given by the equation, 

{p(vxV z ))(h ) = (p t (v x v z ) 1 -){h ) + (p l (v x v z )i)(h ), (34) 

which describes the difference between the horizontal mo- 
mentum of particles impacting onto the sediment bed 
per unit time and soil area (— (pi(vxV z )i)(h )) and the 
horizontal momentum of particles leaving the sediment 
bed per unit time and soil area ((p-\-(v x v z }-f)(h )). This 
momentum change per unit time and area is the conse- 
quence of collisions between particles within the sediment 
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bed (z < h )- Thus, (p(v x v z ))(h ) is an effective fric- 
tional force which the soil applies on the transport layer 
per unit soil area. Bagnold [31, 35] was the first who 
proposed that this frictional force is approximately pro- 
portional to the normal force which the transport layer 
applies on the sediment bed (Coulomb friction law), inde- 
pendently of whether the transport regime is subaqueous 
or aeolian. Models for saturated sediment transport us- 
ing this Coulomb friction law have been successfully vali- 
dated through comparison with experiments, thus giving 
support to the Coulomb friction law adapted to sediment 
transport (e.g. [9, 31, 35, 36]). It seems reasonable that 
the Coulomb friction law should be also approximately 
valid in situations out-of-equilibrium [16], provided the 
sediment flux is close to its saturated value. 

The Coulomb friction law can be used in order to ex- 
press the second term on the right-hand-side of Eq. (28) 
in terms of M and V. Since the buoyancy-reduced grav- 
ity force is typically much larger than the vertical drag 
force, and since the vertical drag force acting on a particle 
in its upward motion (f) tends to cancel the vertical drag 
force acting on the same particle in its downward motion 
(J.) , the normal force which the transport layer applies on 
the sediment bed can be taken as approximately equal to 
the buoyancy-reduced gravity force [9, 31] corrected by 
the added mass force. Hence, the second term on the 
right-hand-side of Eq. (28) can be written as, 



(p(v x v z })(ho) 



gM, 



(35) 



where g = (s — l)g/s is the buoyancy-reduced value of 
the gravity constant, g, and /i is the Coulomb friction 
coefficient. 



3. The conservation equations in their final form 

By substituting Eqs. (31) and (35) into Eq. (28), we ob- 
tain the momentum conservation equation in terms of M 
and V, whereas Eq. (27) gives the mass balance. There- 
fore, the mass and momentum conservation equations in 
their final form read, 



d(MV) 



dx 

d(MV 2 ) _ 3M 



(36) 



d.r 



4c a sd 



C d (V r ) ■ V 2 - —gM. (37) 



We note that Eq. (37) can be used to obtain the saturated 
value V rs of the velocity difference V r . By using d/dx = 
(saturated sediment transport), we obtain, 



— • C d (V rs ) ■ V 2 S = fxg, 
which can be numerically solved for V rs . 



(38) 



III. OBTAINING THE FLUX SATURATION 
LENGTH OF SEDIMENT TRANSPORT 

In this Section, we use the results presented in last 
Section in order to derive a closed expression for the sat- 
uration length as a function of the attributes of sediment 
and flow, both for aeolian and subaqueous regimes. The 
derivation of the saturation length equation is the subject 
of Section III A. In Section III B we present and discuss 
the resulting equation for the saturation length. In Sec- 
tions III C and III D we show how the saturation length 
equation can be applied to compute L s in the aeolian and 
subaqueous regimes, respectively. 



A. Derivation 

Close to equilibrium, M and V saturate simultaneously 
following a certain function M(V), where M s = M(V S ). 
This function is linked to the characteristics of the ero- 
sion and deposition of bed material and thus to the un- 
known shape of (p(v z ))(h ) as a function of M and V in 
Eq. (36). Moreover, also the mean fluid velocity U will 
saturate following a certain function U(V) close to the 
saturated regime, since U is influenced by the feedback 
of the sediment transport on the fluid flow. Therefore, 
Eq. (32) becomes, 



V r (V) = U{V) - V. 



(39) 



By taking into account that both M and U are func- 
tions of V, and by using Eq. (38), we can rewrite the 
momentum balance Eq. (37) in such a way to obtain the 
following expression for dV/dx, 



dV_ 

dx 



il(V) = A(V) ■ B[V r {V)], 



(40) 



where the functions A(V) and B(V) are given by the 
equations, 



A(V) 



3M(V) 



(41) 



4sdc a c v (2V ■ M(V) + V 2 ^§p-y 
B{V r ) = C d {V r ) ■ V r 2 - C d (V rs ) ■ V 2 S . (42) 
Furthermore, since Q(V) = M(V)V, we obtain, 



dV 
dQ 



00-("00 + V^2)-\ (43) 



In this manner, using Eq. (40), r(V) can be written as, 

T(V) = ^(V) = (M(V) + V^1)^V).(U) 

Using Eqs. (43) and (44), we can write Eq. (2) for the 
saturation length as, 



L„ = - 



dT dV 
dVdQ 



v=v s 



dn 

dV 



(45) 



v=v s 
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where wc further used that fl(V s ) = 0. 

Calculating L s through Eq. (45) requires obtaining an 
expression for dCi/dV, where Sl(V) is defined in Eq. (40). 
However, fi(V) incorporates, through the function B(V) 
defined in Eq. (42), a dependence on the equilibrium 
value of the relative velocity V r , i.e. V rs . In order to 
obtain an expression for V rs , we solve Eq. (38) for V rs 
using the drag law of Julicn [37] for natural sediment, 
which writes, 



Cd(v r ) = — - 



1.5, 



(46) 



whereas we find that the specific choice of the drag 
law has a negligible effect on the value of L s obtained 
from our calculations. By substituting the expression for 
Cd(V rs ), obtained with Eq. (46), into Eq. (38), and solv- 
ing this equation for V rs , yields, 



x/8psgd 3 + (24^) 2 - 24u 



3d 



This equation is, then, used to compute B(V r ) through 
Eq. (42), whereupon £l(V) can be obtained using 
Eqs. (39) and (40). The resulting expression for the sat- 
uration length, computed with Eq. (45), reads, 



L„ = 



(2 + c M ) V rs V s F 



(J>9 



where the quantity, 



cm 



V s dAI 
Ms dV 



(Vs), 



(48) 



(49) 



describes the relative change of M with V close to the 
saturated regime, while F is given by, 



F = C d (V rs ) ■ V rs 



fd(C d V r 2 ) 



V dV r 
24V ra u/d + 1.5V? s V T 



Vr=V r , 

+ Vovjd 



2W rs v/d + 3V r 2 s 2V rs + 16i//d' 



(50) 



and thus F encodes information about the drag law. 

In order to obtain our final expression for L s , we need 
to express 4K(V^). We note that, for the saturated state, 
the mean flow velocity U is dominantly a function of the 
shear velocity and the shear velocity at the bed [9, 27], 
that is, 



u b = \jTf(h )/pf. 



(51) 



The shear velocity at the bed, Ub, is reduced due to 
the feedback of the sediment transport on the fluid flow, 
where Tf(z) is the fluid shear stress profile. We can ex- 
press Ub using the inner turbulent boundary layer approx- 
imation of the Navier-Stokes equations. These equations 
approximate the Navier-Stokes equations for heights z 
much smaller than the height 5b of the boundary layer, 
which is the region in which we are interested since typ- 
ically z C ut <C 5b (Eq. (6)). George [38] derived the in- 
ner boundary layer approximation of the Navier-Stokes 



equations in the absence of an external body force. In 
the presence of an external body force, these equations 
must be slightly modified by adding the body force term 
in the momentum equations. The horizontal momentum 
equation thus writes [38], 



dry 
dz 



-F a 



a;body 5 



(52) 



where F x b d y (z) is the horizontal body force per unit 
volume acting on the flow at each height z. i^body is 
the drag force per unit volume which the particles apply 
on the fluid. In other words, i^body is the reaction force 
per unit volume of the horizontal force per unit volume 
which the fluid applies on the particles. That is, 



-^cbody (F- 



(53) 



^ r f (z cut ) 



We then substitute Eq. (53) into Eq. (52) and integrate 
this equation from z = h Q to z = z cut , thereby using 



p f i 



and j:r P (F x ) = r P {F x 



M(F X ). This leads to, 



Tf 



M 



(F x ) 



(54) 



By substituting this equation into Eq. (51) and using 
Eq. (31), we obtain the following equation for u&, 



Fr 



3 M 



u b 



C d {V r )V? 



(55) 



PfK V Pf u * 

Since it* does not depend on V, we can now express 
du , 

dV 1 



if(V s ) as, 



dU 
d~V 



V s + V rs dUb , T . . 

c u — — (Vs), 



Ubs 



dV 



(56) 



where Ub s is the value of Ub in the equilibrium, and the 
quantity c\j is given by the equation, 



Ubs dU 



(Ubs) 



u bs 



dU 



(ubs)- (57) 



U s du b V s + Vrs du b 

where we used, V rs = U s — V s ■ We note that cjj describes 
the relative change of U with Ub close to the saturated 
regime. Moreover, the derivative ^Tr(V s ) can be calcu- 
lated using Eq. (55) with, M = M(V), V r = V r (V), and, 



3A/ S 



C d (V rs ) ■ V 2 



Pf U * 

which follows from Ub(V s 



1 - ^f, (58) 



dub 
dV 



(Vs) 



Ubs 
1 - 



We thus obtain, 



2u 



bs 



$(Vs. 



FV r! 



CM 
Vs 



(59) 



This expression is substituted into Eq. (56), whereas the 

resulting equation is then solved for (l — ^ (Vs)) 
this is the term involving dU/dV which we need to com- 
pute L s in Eq. (48). In this manner, we finally obtain 
a closed expression for the saturation length, which we 
present and discuss in the next subsection. 
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B. The saturation length equation 

The equation for the saturation length, L s , reads, 

_ (2 + c M )c a c v V rs V s FK 
J-'s — : 1 (.oUJ 

where V rs and F are calculated using Eqs. (38) and (50), 
respectively. In addition, the last factor K on the right- 
hand-side of Eq. (60) is given by the equation, 



K 



1 + 


'cu ■ {v s + v rs y 

2FV rs 


'(I" 1 ) 


1 + 


cuc M ■ (V s + V rs 


)" 


■d-0 


2V S 



i + 


'cu ■ (v s + v rs y 

2FV rs 






l c v cm ■ (y s + v rs 


) 





(61) 



and K thus encodes information about the saturation of 
the transport-flow feedback. In fact, if the saturation of 
the transport- flow feedback is neglected (U = U s ), it fol- 
lows cjj = and thus K = 1. In Eq. (61), we assumed 
that the saturated shear velocity at the bed (itf, s ) and the 
bed shear stress in equilibrium (t/ s (/i )) approximately 
equal u t and r t , respectively, i.e. the minimal shear ve- 
locity and the minimal shear stress at which sediment 
transport can be sustained, 



T fs (h ) 

Ubs 



n, 

u t . 



(62) 
(63) 



In the following, we present arguments which justify this 
assumption. 

For acolian sediment transport, Eqs. (62) and (63) 
are known as "Owen's hypothesis" . These equations are 
known to be approximately valid when u* is close to the 
threshold (e.g. Figure 2.10 in Ref. [6]). However, as 
it* increases, u^ s actually decreases away from u t [6, 9]. 
Nonetheless, the approximation which we use in Eq. (61) 
is reasonable even for large shear velocities, since, when 
u* is significantly larger than ut (which means > 2u t 
for Earth conditions with cm — cu = 1), we have that, 



K 



V s 



CmFVts 



(64) 



which is nearly independent of Uf, s . 

For subaqueous sediment transport, Eqs. (62) and (63) 
are known as "Bagnold's hypothesis". This hypothesis 
is widely used in the literature (e.g. [31, 35, 36, 39]), 
although some studies have questioned it (e.g. [25, 40]). 
However, there is evidence from recent studies that this 
hypothesis is approximately fulfilled. In order to review 
this evidence, we use Eqs. (31), (38), and (54) to express 
M s as, 



Ms 



Ca_ 



T fs (h )} 



(65) 



To our knowledge, the only study in which M s has been 
measured as a function of r is the recent study of La- 
jeuncsse et al. [41], who obtained, using video- imaging 
techniques, that, 



M, 



1 



0.4155 



(66) 



Therefore, if we assume Tf s (h ) = r t as in Eq. (62), then, 
by comparing Eqs. (65) and (66) with c a = 1.19 valid 
for subaqueous sediment transport (cf. Eq. (33) with 
s = 2.65), we obtain /j,/c a = 0.415 and thus \i ~ 0.493. 
Indeed, values within the range between /x/c a = 0.3 and 
fi/c a = 0.5 — and thus consistent with the value of /j, es- 
timated above — have been reported from measurements 
of particle trajectories in the subaqueous sediment trans- 
port [42-44]. Further evidence that Bagnold's hypothe- 
sis is approximately correct was provided by the recent 
numerical study of Duran et al. [27]. These authors sim- 
ulated the dynamics of both the sediment bed and of 
transported particles at the single particle scale. Duran 
et al. [27] found that M s g oc (r — it), which is similar 
to Eq. (66) and can satisfactorily explain all simulated 
data with a single proportionality constant. Moreover, 
the authors also found that tj s reduces to r t at a height 
z very close to the top of the bed, z ~ h a (private commu- 
nication with Duran et al. [27]). Given these separate 
lines of evidence, we believe that Bagnold's hypothesis 
is a reasonable approximation. Moreover, we emphasize 
that our analysis for subaqueous sediment transport is 
not affected by this approximation, since we estimate in 
Section III D 2 that cu « and thus K « 1 regardless of 
the value of Ub s - 

In summary, the saturation length of sediment flux, L s , 
can be calculated using Eq. (60), where V rs and F are 
given by Eqs. (47) and (50), respectively, while Eq. (61) 

is used to compute the term, (l — ;nH^s)) > which ap- 
pears on the right-hand-side of Eq. (60). These equations 
include certain quantities which depend on the character- 
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istics of the sediment transport and thus on the transport 
regime. These quantities are c v , cm, cjj, zx, the saturated 
particle velocity V s , and the threshold shear velocity, u t . 
We estimate these quantities for the aeolian regime of 
sediment transport in Section III C and for the subaque- 
ous regime in Section HID. 

C. The saturation length of aeolian sediment 
transport 

In this section, we estimate the parameters c v , cjj , 
Cm, /i, and express the saturated particle velocity V s 
and the threshold shear velocity u t for aeolian sediment 
transport. Note that we estimate these parameters only 
roughly, which is sufficient in the light of the large scatter 
(factor 2 — 4) of the experimental data [18, 19]. 

1. The parameter c v 

The parameter c„ (Eq. (29)) occurs as a prcfactor 
in Eq. (60), and thus determines the magnitude of L s . 
Since, 

<(^-(^)) 2 )=M)-(^) 2 >0, (67) 



oo oo 




we conclude that c„ must be larger than unity, that is, 



c v = > 1. (68) 

K) 

However, experiments on aeolian sediment transport [32] 
show that the change of (v x )(z) with z is small close to h , 
where most of the transport takes place. Consequently, 
the value of c v must be close to unity. 

We estimate c v from experiments on aeolian sediment 
transport [32, 45]. Creysels et al. [32] measured an ex- 
ponentially decaying particle concentration profile, 

p(z) = p(h )e-^ z - h ^l^ (69) 

and a linearly increasing particle velocity profile, 

(v x )(z) = (v x )(h ) + m(z-h ), (70) 

where z p 10mm, (v x )(h ) ~ lm/s, and m w 70s -1 were 
not varying much with it* . Using these measurements, we 
obtain, 



= JJL(l + (V— ^ V^UtM. (71) 

(v x ) 2 \ \(v x ){h )+mz p ) j ( Vx )2 



In order to obtain c v , it remains to estimate (v x )/(v x ) 2 . 
In order to do so, we use measurements of Greeley et al. 
(Fig. 13 of Ref. [45]), who reported a histogram of the 
horizontal particle velocity of the particles located at a 
height Zh = h a + 2cm, from which we obtain, 



v x ) 2 (z h ) 



1.1. 



(72) 



Since the shape of the distribution of the horizontal par- 
ticle velocity does not vary much with the height [46] , we 
thus estimate c„ as, 



1.17= 



(Vx)'' 



1.17- 



Vx) 2 (z h ) 



1.3. 



(73) 



2. The parameter cu 

In contrast to c„, the parameter c\j (Eq. (57)) signifi- 
cantly influences the functional shape of L s as a function 



of u*. cjj characterizes the significance of the transport- 
flow feedback for the saturation of the sediment flux. For 
instance, cu = means that the transport-flow feedback 
does not affect the saturation process since the flow is 
already saturated (U = U s from Eq. (57)). 

In order to estimate cjj, we need to know how the 
mean fluid speed U behaves as a function of the feedback- 
reduced bed shear velocity Ub (see Eq. (57)). For aeolian 
sediment transport, the fluid speed is strongly suppressed 
by the reaction drag forces which the transported grains 
apply on the wind. The feedback is, in fact, so strong 
that the mean fluid speed in the transport layer changes 
only weakly with u* [9] . In leading-order approximation, 
the mean fluid speed is thus proportional to Ub, 

UsUb (74) 



U 



We thus obtain, from Eq. (57), 
Cu ~ 1. 



(75) 



We note that a value of cjj close to unity is obtained even 
if the more complicated dependence of U on Ub, obtained 
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from modeling saturated sediment flux [9], is taken into 
account. Eq. (75) is approximately valid for «* < 4it t [9]. 
Beyond this range, turbulence- induced fluctuations of the 
shear velocity, neglected in the present work, should af- 
fect the value of cjj. 



3. The parameter cm 

The parameter cm, given by Eq. (49)), occurs as a 
prefactor in Eq. (60) and it further affects the functional 
shape of L s (u*) since cm itself affects the feedback term, 
(1 — w(^s)) -1 . cm encodes the relative importance of 
the respective relaxation processes M — > M s and V — > 
V s for the saturation of the sediment flux. There are 
two extreme cases: cm = and cm — > oo. The case 
cm = means that the saturation of M towards M s is 
much faster than the saturation of V towards V s , while 
the opposite situation corresponds to the case cm — > oo. 

In order to estimate Cm for the aeolian regime of sedi- 
ment transport, we first estimate how the function M(V) 
behaves close to the saturated regime. For this purpose, 
we make use of the fact that, for aeolian sediment trans- 
port, the dominant mechanism which brings grains of the 
sediment bed into motion is the ejection of bed grains due 
to impacts of already transported grains, a mechanism 
known as "splash" (see e.g. [9, 47]). It is known that 
ejection of new grains is mainly due to the impacts of 
the fastest transported particles, whereas the impacts of 
slow particles have a negligible effect on the splash pro- 
cess [47, 48]. Indeed, the speed of a fast impacting grain 
mainly determines the number of ejected grains, but not 
their ejection velocities, as found in experiments [49, 50]. 
The ejected particles are typically slow compared to the 
rebound speed of the impacting particle. In other words, 
the impact of a fast grain naturally results in two species 
of particles: a single (fast) rebounding particle and many 
ejected (slow) particles. 

Using numerical simulations of splash and particle tra- 
jectories, Andrcotti [48] could observe these two distinct 
species in the characteristics of transported particles. 
The author noted that the slow species ("reptons") ac- 
counts for the majority of transported mass per unit soil 
area (M). Furthermore, the author's analysis suggested 
that the impact flux of reptons, and thus, in good ap- 
proximation, the total transported mass M, adjusts to 
changes of the impact flux of the fast species ("saltons") 
within a distance much shorter than L s . Therefore, it 
seems reasonable to treat M as locally equilibrated with 
respect to the impact flux of saltons. The locally equi- 
librated value M eq of M is proportional to the number 
of ejected particles per impact, which in turn is propor- 
tional to the impact speed of saltons [47] and thus ap- 
proximately proportional to V. A rough estimate of the 
function M(V) is therefore, 

M«M eq «^, (76) 



which yields, 

cm « 1. (77) 

4- The Coulomb friction coefficient, \i 

The Coulomb friction coefficient, fi, occurs as a pref- 
actor in Eq. (60), and also changes the functional shape 
of L s through Eq. (47). fi can be determined indirectly 
from measurements of the saturated mass of transported 
particles M s as a function of the shear velocity u* , which 
fulfills the equation [9, 27], 

M s = C -^J-.[ul-ul], (78) 

Note that this equation is Eq. (65) with Tf s (h ) = r t 
(Eq. (62)). The value, 

M~l, (79) 

was found in a previous work [9] through determining /i 
indirectly both from experiments as mentioned above and 
from numerical simulations of aeolian sediment transport 
in the equilibrium. 

5. The saturated particle velocity V„ 

The saturated particle velocity V s is dominantly con- 
trolling the dependence of L s on in Eq. (60). Since 
the dependence of V s on is rather weak for aeolian 
sediment transport [9], the saturation length L s will not 
change much with it*. Here we use an expression for 
V s which has been obtained in a recent work [9], since 
the values of saturated sediment flux obtained using this 
equation produced excellent quantitative agreement with 
measurements [9] . The expression for V s reads [9] , 

V. = V t + ^.l*(%) + ±.FJ*), (80) 
where Vt and F 1 {x) are given by the equations, 

Vt = Vo + vV rs w . th Vo = lg 2 L^JK_ (81) 

1 - r] y nppd 

F 7 (x) = (1 - x) ■ ln(1.78 7 ) + 0.5 • (1 - x 2 ) ■ Ei( 7 ) 

+ 1.154- (1 + a;lna;) • (1 - a;) 2 ' 56 . (82) 

In these equations, Ei(a;) is the exponential integral 
function, k — 0.4 is the von Karman constant, and 
( = 5 x 10" N/m is a dimensional parameter encoding 
the influence of cohesion, while rj = 0.1 and 7 = 0.17 are 
empirically determined parameters [9]. The saturated 
particle velocity V s for transport in the aeolian regime 
can be obtained by itcrativcly solving Eq. (80) for V s and 
using the expressions for V t and F^(x) given by Eqs. (81) 
and (82), respectively. 
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6. The threshold shear velocity u t 

We calculate the threshold shear velocity u t by using 
the following equation, which has been obtained from an 
analytical model for aeolian sediment transport in the 
equilibrium [9], 



k ■ (V rs + V ) 

(1 -77) -In^Zmt/Zo)' 



where z m t is given by the following equation [9], 



(83) 



(84) 



In the equation above, (3 = 0.095 is an empirically de- 
termined parameter [9], while z , which is the surface 
roughness of the quiescent sediment bed, is given by the 
equation [9, 51], 

z = dexp(— kB), with, (85) 
B = 8.5 + (2.51ni? p - 3) cxp [-0.11(lni? p ) 2 5 ] , 

where R p = We note that the ratio between u t 

(which is the threshold for sustained transport) and the 
fluid threshold un required to initiate transport in the 
aeolian regime depends strongly on the environmental 
conditions. Eq. (83) yields u t w 0.8uft for Earth con- 
ditions, which is in agreement with measurements [1]. 
However, the ratio Wt/itft under Martian conditions can 
be as small as 10%, as also found from numerical simula- 
tions [6, 52]. Indeed, Eq. (83), which was obtained from 
the same theoretical work leading to Eq. (80), has been 
validated by comparing its prediction with outcomes of 
numerical simulations [47] under a wide range of fluid- 
to-sediment density ratio and particle diameter, thereby 
leading to excellent quantitative agreement (see Fig. 13b 
of Ref. [9]). 



transport under water using particles of average diameter 
d = 2.24mm and relative shear velocity u*/itt = 2.1 [41]. 
In these experiments, particles were considered as being 
transported if they had a velocity larger than a certain 
cut-off value, v c [41]. The distribution of horizontal ve- 
locities for these transported particles was fitted using an 
exponential distribution, 



1 



Pv{Vx) = TF eX P 
V f 



V X ~ Vr 



Vf 



(86) 



where Vf ~ llOmm/s. By using this distribution, we 
can compute c v as, 



/ v%P v (v x )dv x 



J v x P v (v x )dv x 



1+ (I + I7 



(87) 



Lajeunesse et al. did not report specific values of v c 
corresponding to specific measurements [41]. Instead 
they mentioned that v c lies within the range between 
10 mm/s and 30 mm/s, depending on the water flow rate. 
Since d = 2.24mm and u*/u t = 2.1 (which are the 
values reported for the measurement of P v (v x )) corre- 
spond to intermediate values for d and u*/ut investigated 
in the experiments [41], we use the intermediate value 
v c = 20 mm/s as an approximate estimate for the aver- 
age cut-off velocity. Using this estimate for v c , Eq.(87) 
yields, 



c v « 1.7, 

for transport in the subaqueous regime. 



2. The parameter cu 



D. The saturation length of subaqueous sediment 
transport 

In this section, we provide expressions for the param- 
eters c v , cjj and cm, as well as for the Coulomb friction 
coefficient, /1, the saturated particle velocity, V s , and the 
threshold shear velocity, u t , for transport in the subaque- 
ous regime. We remark that we estimate these quantities 
only in a rough manner, consistent with the large scatter 
(factor 2 — 4) of the experimental data. 



In contrast to the aeolian regime, the suppression of 
the fluid flow due to the sediment transport in the sub- 
aqueous sediment transport is weak [27]. The mean fluid 
speed U is thus mainly a function of the shear velocity 
it, and the dependence of U on Ub and thus on V is neg- 
ligible. By neglecting this dependence, we obtain, 

cu w 0, (89) 

which is consequence of Eq. (57) with dU /dub ~ 0. 



1. The parameter c v 

We can estimate c v for transport in the subaqueous 
regime from measurements of the distribution P v (v x ) of 
horizontal velocities v x in subaqueous sediment trans- 
port in the equilibrium. Such measurements were un- 
dertaken by Lajeunesse et al. in experiments of sediment 



3. The parameter cm 

In order to estimate cm for subaqueous sediment trans- 
port, we use evidence provided by the recent numerical 
study of Duran et al. [27]. As mentioned before, these 
authors simulated the dynamics of both the transported 
particles and the sediment bed at the single particle scale. 
Duran et al. [27] found that, during flux saturation in 
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subaqueous sediment transport, M changes within a time 
scale which is more than one order of magnitude larger 
than the time scale in which Q changes. This observation 
can be mathematically expressed as, 



AM 
df 



dQ 



df 



and thus, 



Eq. (91) further implies that, 
V dM 



dM 








< 




dt 


dt 



M dV 



<1, 



and thus 



\c M \ < 1, 



(90) 



(91) 



(92) 



(93) 



where we used the definition of cm, which is given by 
Eq. (49). Hence, we estimate cm as, 



CM ~ 0. 



(94) 



However, we note that our model predictions are consis- 
tent with experiments even if we assume a coupling of M 
to V which is as strong as in the aeolian regime — that 
is, even by assuming cm = 1 and thus increasing L s by 
a factor of 1.5 as compared to the value obtained with 
cm = [53] . This means that the saturation length in the 
subaqueous regime is not very sensitive to the value of 
cm within the range between and 1 (whereas the latter 
value corresponds to sediment transport in the aeolian 
regime) . 



4- The parameter fi 

As obtained in experiments on subaqueous sediment 
transport in the equilibrium, the average mass flux M s 
approximately follows the expression [41], 



0.415,9 



t - r t 



(95) 



By comparing this equation with Eq. (65) with T/ = T t 
and c a = 1.19 (see section IIIB), we obtain, 

H « 0.493. (96) 

for sediment transport in the subaqueous regime. 

5. The saturated particle velocity V s 

It has been verified in a large number of experimen- 
tal studies [31, 35, 40, 41, 54-56], that the equilibrium 
particle velocity in the subaqueous regime of transport 
approximately follows the expression, 




Fig 



Average value of the dimensionless fluid speed J&j 



— Jf£- , as a function of the dimensionless shear velocity . 

sga 7 sga 

For the symbols, the average dimensionless particle speed 
was obtained from measurements [41], while we computed 
V rs using Eq. (47) with fj, = 0.493. The black solid line corre- 
sponds to the best fit to the experimental data using Eq. (97), 
which yields a « 4.6. 



where a is a dimensionless number. We note that the 
above expression is consequence of the equation, V s = 
U s — V rs , where U s is taken proportional to u*. In order 
to obtain V s for sediment transport in the subaqueous 
regime using Eq. (97), we calculate V rs using Eq. (47) 
and use the value a ~ 4.6, which we have obtained by 
comparing the prediction of Eq. (97) with measurements 
of V s as a function of from experiments on subaqueous 
sediment transport in the equilibrium [41] (see Fig. 1). 



6. The threshold shear velocity ut 

The threshold velocity for sustained sediment trans- 
port, ut, in the subaqueous regime is computed by using 
the equation, 



lit 



(98) 



where the threshold Shields parameter O t is obtained 
through an empirical fit to the Shields diagram [57] . The 
resulting expression for Ot reads [57], 



e t = 



0.273 
1 + L2Dl 



0.046- l- 0.576e 



-0.02D, 



) , (99) 



V. 



an* 



(97) where, £>* = d^/s^Ji^. 
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IV. DEPENDENCE OF THE SATURATION 
LENGTH ON PARTICLE SIZE AND FLUID 
SHEAR VELOCITY 

In order to understand the morphodynamics of sedi- 
ment landscapes under water and on planetary surfaces, 
it is important to understand the behaviour of the flux 
saturation length as a function of the relevant attributes 
of sediment and fluid. In particular, the size of planetary 
dunes can serve as a proxy for the saturation length of 
extraterrestrial dune fields, which can be used to infer 
the local fluid shear velocity (u*) and average size (e£) of 
the constituent sediment [5, 6, 58, 59]. In fact, the depen- 
dence of L s on u* has been subject of intense debate in 
previous theoretical works [16, 18-20, 58]. It is therefore 
useful to perform in this Section a systematic study of 
the saturation length as a function of these two relevant 
parameters under different environmental conditions. 



A. The saturation length as a function of the fluid 
shear velocity, m* 

Fig. 2 shows the dependence of L s /(sd) on u„/ut for ae- 
olian sediment transport on Earth (brown solid line) and 
Mars (red dashed line) and for sediment transport under 
water (blue dash-dotted line) computed using Eq. (60) 
for particles with mean diameter d = 250/.tm. The be- 
havior of L s with predicted from Eq. (60) is in clear 
contrast to the scaling relation L s w 2sd proposed in 
previous works [17-20]. This approximate scaling, which 
includes no dependence of L s on it* , was obtained by as- 
suming that the acceleration of transported particles due 
to fluid drag is the dominant relaxation mechanism, and 
by neglecting the entrainment of sediment bed particles 
due to fluid lift as well as the entrainment of sediment bed 
particles and the deceleration of transported particles re- 
sulting from grain-bed collisions. In our more compre- 
hensive model for saturation of sediment flux, however, 
all these aforementioned relaxation processes are taken 
into account. There are two main reasons for the dis- 
parity in the behaviour of L s with observed in our 
model and in the model of Refs. [17-20]. 

First, our expression involves a significant dependence 
of L s on u*/ut due to the dependence of L s on the av- 
erage particle velocity V s and the feedback term K (see 
Eq. (61)), both of which are functions of u*/u t . In partic- 
ular, for the subaqueous regime, V s is a strongly increas- 
ing function of u*/ut thus explaining the strong increase 
of L s with u*/ut in this regime. Furthermore, we see in 
Fig. 2 that in our model the dependence of L s on it, in the 
aeolian regime is small, but not negligible as suggested in 
the model of Refs. [17-20]. Indeed, the dependence of L s 
on and K is a consequence of considering grain-bed 
collisions and the transport-flow feedback, respectively, 
for the saturation of the sediment mass flux Q — both 
neglected in the models of Refs. [17-20]. Second, in con- 
trast to the models proposed in these works, our model 
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Fig. 2. Dimensionless saturation length L 3 /(sd) versus di- 
mensionless shear velocity u,/u t for particles with mean di- 
ameter d = 250/im, computed using Eq. (60). The brown, 
solid line corresponds to aeolian sediment transport on Earth 
(p p = 2650kg/m 3 , p w = 1.2kg/m 3 , g = 9.81m/s 2 , v = 
1.5 x 10 _5 m 2 /s), the red, dashed line corresponds to aeo- 
lian sediment transport on Mars (p p = 3000kg/m 3 , p w = 
0.0185kg/m 3 , g = 3.71m/s 2 , v = 6.4 x 10" 4 m 2 /s), and the 
blue, dash-dotted line corresponds to subaqueous sediment 
transport under water (p p — 2650kg/m 3 , p w — 1000kg/m 3 , 



g = 9.81m/s 2 , 
line shows L 3 



v — 10 m /s). Moreover, the black, dotted 
= 2sd as proposed by Refs. [17-20]. 



considers the dependence of the drag coefficient Cd on the 
particle Reynolds number, R ep = V rs d/v (see Eq. 31). As 
can be seen in Fig. 2, the difference between the parti- 
cle Reynolds numbers on Earth (R ep w 30) and Mars 
(R ep ~ 1) results in an order of magnitude difference 
between L s /(sd) on these two planetary bodies, which 
occurs because the normalized particle velocity V s j 'y/sgd 
increases strongly with R ep when R ep is of order unity 
[9]. 

It can also be seen in Fig. 2 that the saturation length 
L s on Mars decreases with w* when it* is sufficiently close 
to ut, even though V s increases with it, in this regime. 
This surprising behaviour is a consequence of the feed- 
back term K, which, for Mars conditions, decreases with 
sufficiently close to u t and thus overcompensates the 
tendency of L s to increase with V s . In contrast, for Earth 
conditions, the feedback term K increases with it* close 
to Ut- This qualitative difference in the change of K with 
it, between Earth and Mars conditions can be understood 
by noting that, 



V, 



CAlFV r , 



< 1 



(100) 



for Mars conditions with sufficiently small u*/ut, while 
for Earth conditions (and for Mars conditions with suffi- 
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cicntly large u*/ut) the following relation holds, 



V. 



c M FV rs 



> 1. 



(101) 



The physical origin for the difference in the behaviour 
of K with u* mentioned above lies in the mechanics of 
the reduction of the fluid speed due to sediment transport 
(see Section III A). The fluid velocity in the transport 
layer (U) decreases with the average drag force applied 
by the fluid onto the transport layer (^^C d (V r )V 2 ). 
This average drag force in turn is proportional to both 
the mass density M of transported particles and to the 
acceleration term C d {V r )V 2 . Recalling that cm ~ 1 in 
the aeolian regime, we have that, close to saturation, 
both M and V are smaller (if Q < Q s ) or larger (if 
Q > Qs) than their respective saturated values, M s 
and V s . Hence, if M < M s (M > M a ), it follows 
that C d (V r )V r 2 > C d (V rs )V r 2 s (C d (V r )V r 2 < C d (V rs )V r 2 s ), 
which means that the mass density and the acceleration 
term deviate from their saturated values in opposite di- 
rections. The average drag force applied by the fluid 
onto the transport layer C d (V r )V 2 ) thus can be 

both larger (under Mars conditions for sufficiently small 
u*/ut) or smaller (under both terrestrial and Martian 
conditions for sufficiently large u*/ut) than its saturated 
value. Consequently, the fluid velocity in the transport 
layer (£7) can be smaller or larger than its saturated value 
U s , depending whether Eq. (101) or Eq. (100), respec- 
tively, is fulfilled. Further, if U > U s (Eq. (100)), then 
K < 1 and thus the saturation length L s decreases in 
comparison to the situation in which the saturation of 
feedback is neglected {K = 1). In contrast, if U < U s 
(Eq. (101)), K > 1 and thus the saturation length L s in- 
creases in comparison to K = 1. The deviation of K from 
K = 1 becomes stronger with increasing u*/ut because 
the effect of the sediment transport on the fluid veloc- 
ity increases with it*/u t . This explains why the feedback 
term K can both increase or decrease with u*/ut- 



B. Dependence of the saturation length on the 
average particle size (d) 

Fig. 3 shows the dependence oiL s /(sd) on d for aeolian 
sediment transport on Earth (brown solid line) and Mars 
(red dashed line), and for subaqueous sediment trans- 
port (blue dash-dotted line), computed using Eq. (60) 
for different values of the fluid shear velocity, namely, 
u* = Ut, u* = 2u t , and = 4it t . It can be seen that 
the rescaled saturation length L s /{sd) displays a com- 
plex behaviour with d. The dependence of L a /(sd) on 
d is controlled by two main factors: the dependence of 
the drag coefficient and of the Shields parameter on the 
particle Reynolds number (R ep ), and the dependence of 
the saturation of the sediment flux on the transport-flow 
feedback. The significance of each one of these factors for 
the dependence of the saturation length on the average 



grain diameter depends on the transport regime, as we 
will discuss in the next paragraphs. 

Subaqueous regime of transport — The change in the 
transport-flow feedback due to the saturation of the 
flux is negligible in the subaqueous regime, since in this 
regime cjj and thus K ~ 1. This behaviour explains 
why all curves L s / {sd) versus d in Figs. 3a— c correspond- 
ing to transport under water display the same qualitative 
behaviour independently of u*. In each of these curves, 
L s /(sd) first increases with d, then reaches a maximum, 
and finally approaches a constant value. The origin of 
this behavior is that L s is proportional to V rs V s F (cf. 
Eq. (60)), and thus its dependence on the grain size is 
determined by the behaviour of V rs V s F with d. For suf- 
ficiently large particle diameters d (large R ep ), both the 
drag coefficient C d and the Shields parameter Q t are ap- 
proximately independent of d, and thus both V rs / 'y 'sgd 
and V s /y/sgd are roughly independent of d (see Eq. (47), 
as well as Eqs. (97) and (98) for constant u*/ut), while 
F w 0.5 (see Eq. (50)). Hence, L s /(sd) is nearly in- 
dependent of d for sufficiently large particle diameters 
in the subaqueous regime. However, for smaller parti- 
cle diameters d (smaller R ep ), both C d and 0t, and thus 
V rs I W sgd, Vs/^/sgd and F, incorporate a dependence on 
d. In this regime, V rs / yfsgd increases with d, whereas 
both Vg/y/sgd and F decrease with d. The increase 
of Vrs/y/sgd with d thereby dominates the behaviour of 
L s /(sd) for small particle diameters, while the decrease 
of V s /\/sgd and F with d dominates the behaviour of the 
saturation length with d for large particle diameters. 

Aeolian regime under terrestrial conditions — The 
same qualitative behavior of L s /{sd) with d observed in 
the subaqueous regime occurs for aeolian sediment trans- 
port on Earth at shear velocities close to the thresh- 
old, as can be seen in Fig. 3a. Two main factors dic- 
tate the observed dependence of L 8 /(sd) on d for aeolian 
transport under terrestrial conditions. First, close to the 
threshold, the fluid velocity is almost undisturbed by the 
particle transport, since only a few particles are trans- 
ported. Hence, the feedback term K associated with the 
aeolian regime of transport under terrestrial conditions 
when is close to threshold for sustained transport, lit, 
is close to unity, as it is in the subaqueous regime. Sec- 
ond, V rs I V 'sgd, V s / V sgd and F behave qualitatively in 
the same manner with d as they do in the subaqueous 
regime. Therefore, also L s /(sd) for aeolian transport un- 
der terrestrial conditions with u* close to u t depends on d 
in the same manner as it does in the subaqueous regime. 
In contrast, for large shear velocities, the qualitative be- 
havior of L s /(sd) with d observed for aeolian sediment 
transport on Earth is qualitatively different from the one 
observed in the subaqueous regime (see Figs. 3b— c). This 
is because the saturation of the transport-flow feedback 
for aeolian sediment transport plays a relevant role for 
large u*. In this regime, the approximation K f=a 1 is 
not valid anymore, instead the feedback term K follows 
Eq. (64). Consequently, L s is proportional to V 2 and 
not to V rs V s F. Since V s / y/sgd decreases with d before it 
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Fig. 3. Normalized saturation length L s /(sd) versus mean particle diameter d for three shear velocities, namely, u* = ut, 
u* = 2u t , and u* = 4it t , computed using Eq. (60). The brown, solid line corresponds to aeolian sediment transport on Earth 
(p p = 2650kg/m 3 , p w = 1.2kg/m 3 , g = 9.81m/s 2 , v — 1.5 x 10 5 m 2 /s), the red, dashed line corresponds to aeolian sediment 
transport on Mars (p p = 3000kg/m 3 , p w = 0.0185kg/m 3 , g = 3.71m/s 2 , v = 6.4 x 10 _4 m 2 /s), and the blue, dash-dotted line 
corresponds to subaqueous sediment transport under water (p p = 2650kg/m 3 , p w = 1000kg/m 3 , g = 9.81m/s 2 , v = 10 _6 m 2 /s). 
Moreover, the black, dotted line shows L s = 2sd as proposed by Refs. [17-20]. 



reaches an approximately constant value, so does L s /(sd) 
as can be seen in Figs. 3b— c. 

Aeolian regime under Martian conditions — For aeo- 
lian transport under Martian conditions, L s /(sd) shows 
a qualitative behavior with d that is different from the 
one of aeolian transport under terrestrial conditions. The 
origin of this discrepancy is a different qualitative behav- 
ior of the normalized average particle velocity V s /y/sgd, 
as we explain in the following. In aeolian sediment trans- 
port, the average velocity of particles impacting and leav- 
ing the sediment bed (V a ) is nearly independent of prop- 
erties of the fluid [6, 7, 9, 32, 46-48, 52, 60]. Rather, 
V is largely controlled by the characteristics of the sed- 
iment bed, for instance by cohesive interparticlc forces. 
These forces increase in importance with decreasing par- 
ticle size [6]. Due to these forces, V /\fgd increases with 
decreasing d (see Eq. (81)). At the same time, as already 
mentioned, V rs /\/sgd increases with d. Both V /\fsgd 
and V rs I V 'sgd control the normalized average particle ve- 
locity V s /y/sgd (see Eqs. (80) and (81)). 

The qualitative difference between aeolian sediment 
transport on Earth and Mars is due to V /V rs being 
approximately 2.5 times larger on Earth than on Mars. 
This difference in the scaling of V /V rs implies that the 
decreasing trend of V j y 'sgd with d has a smaller effect 
on the saturation length for Martian conditions than it 
does for L s under terrestrial conditions. Indeed, while 
the value of V s j 'y 'sgd of aeolian transport under terres- 



trial conditions decreases with d (due to the decrease of 
V 1 y 'sgd with d) , the Martian value of V s /\/sgd increases 
with d. This behaviour, which is opposite to the one of 
terrestrial transport, is due to the increase of V rs /y/sgd 
with d — except for small particle diameters for which 
V rs I V sgd becomes much smaller than V / \fsgd. 

The complex behavior of the saturation length L s plot- 
ted in Figs. 2 and 3 suggests that both the entrainment 
of sediment bed particles by fluid lift and grain-bed col- 
lisions, and the momentum change of saltating particles 
due to drag and grain-bed collisions have a considerable 
influence on the saturation of the sediment flux. Thus, 
these relaxation mechanisms of the sediment flux proba- 
bly cannot be neglected as done in previous studies [16- 
20, 22]. Moreover, the relaxation of U, which has also 
been neglected in previous studies, seems to play a sig- 
nificant role in the saturation length of aeolian sediment 
transport, for which the feedback term K (Eq. 61) does 
not vanish since c\j ~ cm ~ 1- In fact, for sufficiently 
large it*, the feedback term K follows Eq. (64), which 
means that L s scales with V s 2 instead of scaling with 
V rs V s F for small it*. Hence, our study suggests that the 
saturation of U cannot be neglected for aeolian sediment 
transport. 

This scaling of L s with V s 2 for large u* (it* > 2ut for 
terrestrial conditions) has an interesting implication for 
the behaviour of L s with the density of the transported 
particles in aeolian sediment transport. For environmen- 
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tal conditions in which V s is strongly controlled by the 
particle slip velocity V a (for instance terrestrial condi- 
tions), L s depends only weakly on the particle density p p 
since V is independent of p p . This means that the respec- 
tive values of the saturation length rescaled by the grain 
diameter, L s /d, of sand and snow particles in aeolian 
transport under terrestrial conditions are of the same or- 
der of magnitude, even though the density of snow parti- 
cles is almost an order of magnitude smaller than the den- 
sity of sand particles. This surprising prediction, which is 
in clear contrast to the prediction L s w 2sd = 2(p p /pf)d 
of Refs. [17-20], can be checked against estimates of L s 
from the size of snow and sand dunes on Earth. Such 
measurements, indeed, indicate that L s /d is of the same 
order of magnitude for the transport of snow and sand 
particles [53]. 

V. CONCLUSIONS 

In conclusion, we have presented a model for flux sat- 
uration in sediment transport which, for the first time, 
accounts for both relevant relaxation processes of sed- 
iment flux identified in previous works — namely, the 
saturation of the mass density of transported particles 
and the relaxation of particle velocities — as well as for 
the different types of sediment entrainment prevailing un- 
der different environmental conditions. Furthermore, our 
model accounts for the saturation transient of the fluid 
velocity within the transport layer, which is associated 
with the saturation of the transport-flow feedback inher- 
ent to the interaction between the fluid and the particles 
in transport. The main outcome of our analytical treat- 
ment is a closed-expression for the saturation length of 
sediment transport, L s (Eq. (60)), which can be used 
to calculate L s under different environmental conditions 
corresponding to both subaqueous and aeolian regimes 
of sediment transport. In particular, L s predicted from 
our equation is a complex function of the grain diameter, 
d, and of the fluid shear velocity, u*. This behavior is 
in contrast with the scaling of L s with sd [17-20], which 
was obtained from a simplified model that considers only 
the relaxation of particle velocity and thus neglects the 



dependence of L s on u» observed in experiments [21]. 

While the purpose of the present work was to introduce 
our theoretical model for flux saturation and to present 
the analytical derivation of our universal equation for the 
saturation length, in a separate work [53] we show that 
this equation consistently predicts the saturation length 
in different physical environments. Indeed, our satura- 
tion length equation is in good quantitative agreement 
with direct measurements of L s in a wind tunnel, as well 
as with indirect estimates of L s from the size of subaque- 
ous ripples, snow dunes and dunes on Earth, Mars and 
Venus. 

In future studies, our equation can be used to predict 
the scale of dunes under different extraterrestrial envi- 
ronments [53] or to infer attributes of sediment and flow 
in planetary dune fields from the minimal size of barchan 
dunes or from the wavelength of "elementary" dunes 
emerging on dense sand beds. Moreover, our equation 
can be used to calculate the saturation length in aeolian 
dune models [5] and to improve morphodynamic models 
in hydraulic engineering applications in which the satura- 
tion length is usually treated as an adjustable parameter 
[61, 62] . The importance of the saturation length for such 
models has been debated among engineers [63-65], and 
it was concluded that even at fluvial scale the influence 
of sediment transport saturation was significant [65] . We 
thus anticipate that the present work can provide sub- 
stantial contributions to several areas of the geological, 
planetary and engineering sciences. 
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